Infinitesimal propagation equation for atmospheric decoherence with multiphoton 

correlations 
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Previously a set of coupled first order differential equations were derived for the decoherence of 
a pair of spatial mode entangled photons, propagating along different paths through turbulence. 
Here we extend this analysis to the situation where both photons travel along the same path, which 
introduces the effect of multiple photon correlations. The resulting equation now contains additional 
terms that take these multiphoton correlations into account. At the same time, we provide a more 
thorough formulation of the quantized field, starting from a Lorentz invariant formulation, which is 
then explicitly broken by the choice of a particular propagation direction. The effect of the latter 
improvement in the quantization on the form of the final equation is minimal. 

Keywords: Infinitesimal propagation equation, atmospheric scintillation, orbital angular momentum entan- 
glement, decoherence, multiphoton correlation, Lorentz convariant quantization 



I. INTRODUCTION 

While the orbital angular momentum (0AM) states 
of a photon allow higher dimensional representations of 
quantum information flHsj, an entangled photonic state 
that is encoded in terms of 0AM states will suffer degra- 
dation of the encoded information when it propagates 
through a turbulent medium as a result of the decoher- 
ence of the entanglement that is caused by the scintilla- 
tion. To design a quantum communication system one 
therefore needs to understand the decoherence process of 
spatial mode entanglement in a turbulent atmosphere. 

Recently a theoretical framework, in the form of an in- 
finitesimal propagation equation (IPE), has been devel- 
oped to handle this problem Although it extended 
previous work [sI-IMIj which modeled the turbulence with 
a single phase function, the framework in Jl considered 
only single photons going through a particular turbulent 
atmosphere at a time. The effect of multiple photons 
going through the same medium has therefore been ne- 
glected. Here we extend the framework to include the 
additional correlations that are introduced when multi- 
ple photons propagate through the same medium. As a 
result we obtain additional terms in the final expression 
for the IPE. These terms serve to link different photons 
with each other. 

At the same time we also provide a more thorough 
discussion of the quantization of the electric field. Pre- 
viously the quantization of the electric field was ex- 
pressed in the popular way @. However, it can be shown 
that the quantum states that are defined in this way are 
manifestly not Lorentz covariant [lo| (under a Lorentz 
boost the orthogonality relation for these quantum states 
transforms into an expression that is different from the 
original expression). 

On the other hand, one can argue that the final expres- 
sion would inevitable be non-covariant as a result of the 
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explicit violation of Lorentz covariant when one either 
fixes the propagation direction or enforces the field to be 
monochromatic. The question then is whether the final 
expression could be derived consistently from an initial 
Lorentz covariant formulation. In other words, does the 
final expression look the same when one starts with the 
usual non-covariant quantization as it would when one 
starts from a Lorentz covariant quantization and then 
goes, in a consistent manner, through the steps that 
explicitly breaks the Lorentz covariance? We find that 
there is a minor difference between these expressions (the 
integrals in the latter case contain an additional factor of 
1/kz), but in the paraxial limit this difference becomes 
insignificant. 



II. FIELD QUANTIZATION 

In general, quantum basis states are defined in the 
Fourier domain where they represent a plane wave basis 
and depend on the three components of the propagation 
vector k^, ky and kz and on the angular frequency a;. Due 
to the vacuum dispersion relation w = c|k|, any three of 
the four quantities w, k^, ky and kz will fix the fourth. 
Therefore, the phase space integrals for propagating (on 
shell) fields only run over three of these quantities. It is 
natural to select k^, ky and kz as the three integration 
variables, which then fixes the value of w. In this way one 
can define all the quantities in an explicitly Lorentz co- 
variant manner. We'll start with the expressions in this 
form. 

In the current scenario where we consider the decoher- 
ence of 0AM entangled states in atmospheric turbulence, 
we find that the system is not one that evolves in time, 
but rather in space . The quantum states represent ex- 
citations of an optical beam propagating in a particular 
direction through space, which we define as the z-axis. 
In this situation it is more convenient to integrate over 
a;, kx and ky and thereby fix kz. The reason is that the 
system's evolution is considered as a function of the prop- 
agation distance and not as a function of time. Moreover, 
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the input is specified as a function on a two-dimensional 
plane perpendicular to the propagation direction for all 
time (as opposed to a three-dimensional function through 
all space at a fix initial point in time) . Such an input can 
also be expressed in the Fourier domain as a function of 
kx, ky and lu. In the monochromatic case one can fix uj 
and then end up with a two-dimensional function of k^ 
and ky. 

It should be noted that by fixing the z-axis as the prop- 
agation direction, one explicitly breaks the rotation sym- 
metry, which forms part of the Lorentz symmetry. As a 
result our final expression implies an explicit breaking of 
Lorentz invariance. Nevertheless, our initial formulation 
should still be Lorentz invariant to ensure consistency. 
We therefore start with expressions that are defined in 
terms of Lorentz covariant integrals over kx , ky and kz . 
Performing a change of variables, we then end up with 
quantities that are defined in terms of integrals over uj, 
and ky. 



A. Three-dimensional momentum states 

The creation and annihilation operators are defined in 
terms of how they create or annihilate the momentum 
states. The Lorentz covariant orthogonality condition for 
the three-dimensional momentum basis is given by (loj ^ 



(ki,r|k2,s) =w (27r)3j,,, J(ki -k2 



(1) 



where k and the r and s represent 

the spin state, which we include here for the sake of com- 
pleteness. These basis states are generated or destroyed 
by creation and annihilation operators according to 



(k,s| = (0|a,(k) 
|k,s) = V^at(k)lO), 



(2) 



so that they obey the commutation relation 

[a,(ki), 4(k2)] = {2Tz)H,^r ^(ki - k2). (3) 

We define an identity operator that is resolved in terms 
of the three-dimensional momentum basis [loj 



1 = ^1 |k,s)(k. 



d^fc 
UJ (27r)3 



(4) 



^ For a boost along, say, the z-direction, given by k'^ = 7(fcz + 
Pu>/c) and u)' = 7(0; + cf}kz), one finds that the Dirac deha 
function for the z-component transforms as 



dA4 
dfc^ 



where we made use of the vacuum dispersion relation. To avoid 
this non-covariant transformation, the expression for the orthog- 
onality condition needs an extra factor of u). 



where uj is given by 

uj^c{kl + kl + klY'\ (5) 

Integral signs without explicit integration boundaries al- 
ways imply integration from —00 to 00. 

The identity operator defined in Eq. ^ can be used 
to find the expansions of arbitrary one-photon states in 
terms of the three-dimensional momentum states 



Y 
d^fc 



(6) 
(7) 



where the momentum space wave function is given by 

vE'(k,s) = (k,s|vI') (8) 

and * represents the complex conjugate. 

The one-photon states are normalized so that (^'|^') — 
1, which then implies that 



1 = 5^1 \^{\.,s)Y 



UJ (27r) 



(9) 



B. Choosing a propagation direction 

Now we fix the propagation direction to be the z- 
direction and redefine the quantities in terms of the two- 
dimensional momentum-energy states |K,a;,s), instead 
of the three-dimensional momentum states |k, s), where 
K = kxX + kyij and s denotes the spin state. By fixing 
a specific direction for propagation we explicitly break 
rotation invariance and, by implication, also Lorentz in- 
variance. 

From the vacuum dispersion relation cj^ = (?\k\ it fol- 
lows that UJ duj = c^fc, dfcz or 



dfc^ 



Auj. 



(10) 



One needs to be careful with the negative frequencies 
and negative k^s and how they map to each other. In 
the applications that we consider, the angular spectrum 
of the beam only contains nonzero components in the 
positive fc^-direction. Fortunately, since both sides of the 
w-axis map into the positive side of the fcz-axis, we do 
not encounter a problem. 

Applying the change of integration variables given in 
Eq. (Iini) to Eqs. (O and ©, we obtain 



I*) = E / |K,a;,s)vl/(K,a;,s) 

s •' 

(*| =E/**(K,-,.)(K,.,.| (12 



(PK duj 
c^kz (27r)3 

d^K duj 
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where 



1/2 



(13) 



and the momentum space wave function is given by 

^{K,iu,s)^{K,u;,s\^). (14) 

The inverse two-dimensional Fourier transform of the mo- 
mentum space wave function gives the position space 
wave function on the transverse plane at z = 0. 

Using Eqs. (fTTj) . (fT2|) and the orthogonality condition 

(Ki,c.i,r|K2,c.2,s) = c^k, {2nf6r,sS{K,-K2) 

x6{uji-uj2), (15) 

one obtains the normahzation condition of the one- 
photon states in terms of the two-dimensional momen- 
tum states 



i = j2 1 i*(K,c.,s)r 



c^k, (27r)= 



(16) 



C. Monochromatic assumption 

The current application assumes that the optical field 
is monochromatic. For this purpose the momentum space 
wave function is assumed to be given by 



^'(K, s) = c G{K)H{uj - Luo; Slo). 



(17) 



where ojq and Suj are, respectively, the center frequency 
and the (small) bandwidth of the optical field. From now 
on we also drop the spin s, because we assume a uniform 
polarization. To satisfy the normalization requirement 
for the momentum space wave function given in Eq. ([T| 
we define 



H{uj; Suj) — \j exp 



so that 



(18) 



(19) 



This definition of H{uj; Suj) does not actually require the 
field to be monochromatic, unless we take 5uj — >■ 0. If we 
assume that the bandwidth Suj is very small, we can sub- 
stitute u = ujQ inside kz- As a result the normalization 
then reduces to 



where 



1 = I |G(K)f 



kz — f ^ 



d^k 



1/2 



(20) 



(21) 



The monochromatic two-dimensional 
states can be defined by 



momentum 



duj 
2^ 
do; 
2^' 



|K) = J |K,cj)if(w-wo;fc) 

(K| = J H{uj-ujo;5uj){K,uj\ 

and they obey the orthogonality relations 

(K1IK2) =4^'fc. (5(Ki-K2) 

Using Eqs. (f22|) and ([23l) we can define the monochro- 
matic one-photon states as follows 



(22) 
(23) 

(24) 



vf) = 1 |K)G(K) 
G*(K)(K| 



d^K 

Air'^kz 
d^K 
An'^kz 



(25) 



(26) 



in terms of the monochromatic two-dimensional momen- 
tum states. Note that we inserted c's into the definitions 
of the monochromatic two-dimensional momentum states 
in Eqs. ([22]) and ((231) and into Eq. ^7} so that the final 
definitions of the one-photon states in Eqs. ([25]) and (pS)) 
are without c's. 



III. DENSITY OPERATOR IN THE OAM BASIS 

The set of all Laguerre-Gaussian (LG) modes forms 
a complete orthonormal basis. These modes are distin- 
guished by two indices: an azimuthal index i and a radial 
index r. For notational convenience these two indices are 
here combined into one index m = {(,r}. Unless stated 
otherwise, each index used in the subsequent derivation 
always represents both the indices associated with a par- 
ticular LG mode. 

The LG modes are eigenstates of the rotations on the 
transverse plane. Since rotation invariance represents 
conservation of orbital angular momentum (OAM), these 
LG modes are also OAM eigenstates. In fact, the amount 
of OAM of an LG mode is proportional to the azimuthal 
index £ of that mode. The LG modes therefore form an 
OAM basis. 

For the purpose of this derivation we first consider a 
single photon and then generalize the result for the case 
of two photons. The density operator of an arbitrary 
single photon state can be expressed in the OAM basis 
by 



P - 



E 



to) p„ 



(27) 



Each of these OAM states can be expanded in terms 
of the monochromatic two-dimensional momentum ba- 
sis, using Eqs. (gl]) and (pS)) . 



|K)G„(K) 
G;.(K)(K| 



d^K 
in'^kz 

47r2fc, 



(28) 



(29) 
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leading to the following expression for the density op- 
erator in terms of the monochromatic two-dimensional 
momentum basis 



J |Ki) G™(Ki,z)p™_„ 



E 

m,n 

xG:(K„.)(K.|^-^, (30) 

where the dependence on z is shown explicitly to make 
it apparent that this expression for p is only valid on a 
transverse plane for a specific value of z. Here, and also 
later, we use only one integral sign to represent several 
/sT-space integrals. 

By evaluating the summations in Eq. pop . one obtains 
the definition for the density operator in terms of the 
two-dimensional momentum basis. 



p(z) = J |Ki)p(Ki,K2,^)(K2 



, (31) 



where 



p(Ki, K2, ^) - G™(Ki, z)p„,„G;(K2, z). (32) 

The 0AM eigenstates are orthonormal, (m|n) — Sm,n, 
which, together with Eq. (j24p . implies that the momen- 
tum space wave functions G',„(K, z) are also orthonormal 



J g,^(k,z)g;(k,z) 



4:71^ k. 



— , n • 



(33) 



Using Eqs. ((32|) and ([SS)) one can show that 



J g;,(Ki,z)p(Ki,K2,z) 



xG„(K2,2;) 



47r2fc^i 47r2fc;j2 ' 



(34) 



Since the two-dimensional momentum basis and the 
0AM basis are completely equivalent, the definitions in 
Eqs. (|27| and (|3T|) are also completely equivalent and 
Eqs. p2|) and p4|) indicate how one can transform from 
one to the other. 

For two photons the density operator in Eq. ((5D|) can 
be generalized to become 



P(^) 



y |Ki)|K3)G™(Ki,z)Gp(K3,.2) 



E 

m.n,p,q 

xp™,„,p,g g:(K2, z)g;(K4, z)(K2|(K4| 

d^Ki d^K2 d^i^g d^Ki 



47r2fc2i 47r2fc^2 47r2fc^3 47r2fc^4' 



(35) 



For notational convenience we represent the product of 
momentum space wave functions that appear in Eq. (1351) , 
as a single function, 

G„(Ki, z)G:(K2, z)Gp(K3, z)G;(K4, z) 
= F„,„,p,,(Ki,K2,K3,K4,z). (36) 



The expression for p is now given by 



P(^) 



^ ^ Pm,,n.p.q 
m.n,p,q 



|Ki)|K3) 



m,n,p,q 



(Ki,K2,K3,K4,z)(K2|(K4 



d^Ki d^K2 d^K^ d^Ki 
'i'K^kzi 47r2/c22 4:n^kz3 'in'^kz. 



(37) 



Note that inside the expression in Eq. (|37)) the z- 
dependence is carried by the momentum space wave func- 
tions and not by the density matrix elements. This is 
because the transformation of the density operator dur- 
ing propagation over an infinitesimal distance through 
turbulence is caused by the distortion of the momentum 
space wave functions. After such an infinitesimal propa- 
gation these momentum space wave functions no longer 
represent the Fourier transforms of the original modes. 
One needs to re-expand these distorted wave functions 
in terms of the momentum space wave functions of the 
0AM modes and incorporate the expansion coefficients 
in the density matrix elements. Thereby one can transfer 
the z-dependence to the density matrix elements. 

To obtain the full three-dimensional expression for the 
density operator in free-space (without turbulence) one 
can use Fresnel diffraction theory to determine the ex- 
pression at any other value of z. In the presence of tur- 
bulence the expression for the density operator is only 
valid on a specific transverse plane, and it needs to be 
transformed from plane to plane, according to the dy- 
namics of the medium [see Eq. pij) below] . 

Since Fm,n,p,q carries the only z-dependence in the ex- 
pression for the density operator, we focus on how it 
transforms during infinitesimal propagation. At the end 
we apply the transformation to the expression for p. 



IV. EQUATION OF MOTION AND 
INFINITESIMAL TRANSFORMATION 

For a classical electromagnetic field propagating 
through a source free region, the equation of motion, 
which follows directly from Maxwell's equations, is given 
by the Helmholtz equation. 



\/'^E{y.)+n^k^E{y.) = 0, 



(38) 



where n in the refractive index, k is the wave number 
and 'X. — XX -\' yy zz. It is assumed that the polariza- 
tion is uniform and can be ignored, therefore, only the 
scalar part of the electric field i?(x) is considered here. 
The inhomogeneous medium is represented by a spatially 
varying index of refraction 



n = 1 + (5n(x). 



(39) 



This variation is very small {8n ^1), which implies that 
one can approximate the Helmholtz equation as 

V2£;(x) + fc2£;(x) + 25n{x)k^E{x) = 0. (40) 
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Furthermore, we assume that the beam is paraxial and 
propagates in the z-direction. So we define 



E(x) = g(x) exp(— ifcz), 



(41) 



which then leads to the paraxial wave equation with an 
extra inhomogeneous medium term 

V^5(x) - i2kd^g{x) + 25n{yi)k'^ g{x) ^ 0, (42) 

where is the transverse part of the gradient operator. 
The two-dimensional inverse Fourier transform, 



5(x) = J G(K,z)exp(-iK-x) 



4^' 



(43) 



which contains the angular spectrum of the optical field 
G(K, z), is substituted into Eq. (|42]) . to obtain 

d^G(K, z)^^ |K|2g(K, z) - ikN{K, z) * G(K, z), 

(44) 

where ★ indicates convolution and A^(K, z) is the 
two-dimensional Fourier transform of Sn{x.). The z- 
dependence is shown explicitly to indicate that G(K, z) 
depends on the location of the transverse plane along z. 

Note that the (inverse) Fourier transform in Eq. P5|) is 
a purely formal operation on the two-dimensional func- 
tion and therefore does not contain a factor of l/fc^. The 
convolution integral, which comes from the Fourier trans- 
form of the product of two functions, does not contain a 
factor of 1/kz either. 



OAM 
entanglement 

source ,,,,,,,, _ 

Detector 
Detector 

Turbulent 
atmosphere 

FIG. 1: Diagram of an OAM entangled biphoton propagating 
through a turbulent medium toward two detectors. 

The quantum wave function in an interaction-free sys- 
tem obeys the same equation of motion as the classi- 
cal field. As a result, one can identify the momentum 
space quantum wave function with the angular spectrum 
in Eq. (|33]). Thus the expression in Eq. (|44p represents 
the infinitesimal transformation of the momentum space 
wave function during propagation through a turbulent 
(random) medium. It forms the basis of the derivation 
of the IPE. 




V. DERIVATION OF THE IPE 

Here we consider the scenario where two photons (a 
biphoton) , which are entangled in terms of the OAM ba- 
sis, both propagate through a turbulent atmosphere, as 
shown in Fig. [TJ The derivation of the IPE follows the 
same procedure that was followed in The only dif- 
ferences are that the integration measure for the phase 
space integrals contain an extra 1 / fc^-factor and here we 
allow the two photons to go through the same turbulent 
medium, which means that there are additional correla- 
tions that were excluded in 



A. Infinitesimal transformation 

First we consider what happens to Fm.n.p,q after an 
infinitesimal propagation, by setting z z + dz. Ex- 
panding the result to sub-leading order in dz, applying 
Eq. (|44p and taking the limit dz — > 0, we obtain a differ- 
ential equation given by 



dzFm,n,p,q{'i^l, K2, K3, K4, z) 

= ^(|Ki|2-|K2p-H|K3p-|K4n 

X^m,n,p,(j(Ki, K2, K3, K4, z) 

-ik J N(Ki - K, z)F„,„ 



(K,K2,K3,K4,z) 

-iV-(K2 - K, z)i'„,,„^p,g(Ki, K, K3, K4, z) 
+iV(K3 - K, z)i^„.„,p,,(Ki, K2, K, K4, z) 
-N*{K4 - K, z)i^„,„,p,,(Ki, K2, K3, K, z) 

The integral in Eq. (HSl) does not contain the additional 
1/fcz-factor because it is a convolution integral. 

Since the spectra N(K, z) are random functions, 
Fm,n,p,q vly Eq. (|45|) is also a random function. To calcu- 
late the expectation value for F„i.n,p,q we need to employ 
ensemble averaging. (In a slight abuse of notation we'll 
keep on denoting the ensemble average of F, 



m,n^p^q 



smi- 



ply as Fra,n,p,q-) For this purpose we start by integrating 
Eq. (|45p from zg to z, which gives Fm,n,p,q in terms of 
previous versions of itself. Using repeated back substi- 
tution, one obtains a Born series (or Dyson expansion). 
Since each A^(K, z) effectively comes with a factor of G„ 
(the square root of the refractive index structure constant 
G^), higher orders in A^(K, z) are suppressed. However, 
since (A^(K, z)) = 0, one needs to expand the series at 
least up to second order in A^(K, z). The resulting ex- 
pression, without the (A^(K, z))-terms, is given by 



6 



m,n,p,q 



;Ki,K2,K3,K4,z) = F„_„_p^5(Ki,K2,K3,K4,zo) 



'2k 

rzi 



+(-^ - ■^o)7^i:(|Ki|' - |K2|2 + IK3I' - |K4|2)F„,„,p^,(Ki, K2, K3, K4, zo) 



f [ ' /"(A^(Ki - K, zi)iV(K - Ko, z2))F™,„,p,g(Ko, K2, K3, K4, zq) 

J Zn J Zn J 



-(7V(Ki - K, zi)iV*(K2 - Ko, z2))i^™,„,p,g(K, Kq, K3, K4, zq) 

+ (7V(Ki - K, Zi)iV(K3 - Ko, Z2))i^™,„,p,g(K, K2, Kq, K4, zo) 

-(7V(Ki - K, zi)iV*(K4 - Ko, z2))i^,„,„,p,,(K, K2, K3, Kq, zo) 
^(7V*(K2 - K, zi)iV(Ki - Ko, z2))i^,„,„,p,,(Ko, K, K3, K4, zo) 
+ (7V*(K2 - K, zi)iV*(K - Ko, Z2))F™,„.p,,(Ki, Ko, K3, K4, zo) 
^(7V*(K2 - K, zi)iV(K3 - Ko, Z2))F„.„,p,,(Ki, K, Ko, K4, zo) 

+ (7V*(K2 - K, Zi)iV*(K4 - Ko, Z2))i^m,„,p.,(Ki, K, K3, Kq, zo) 

+ (7V(K3 - K, zi)iV(Ki - Ko, z2))^^„.,„,p^,(Ko, K2, K, K4, zo) 
-(iV(K3 - K, zi)iV*(K2 - Ko, z2))^^™,„,p,,(Ki, Ko, K, K4, zo) 
+ (7V(K3 - K, zi)iV(K - Ko, z2))F™,„,p,,(Ki, K2, Ko, K4, zo) 
-(7V(K3 - K, zi)iV*(K4 - Ko, z2))^^™,„,p,,(Ki, K2, K, Kq, zo) 
-(7V*(K4 - K, zi)iV(Ki - Ko, z2))F„,„,p,,(Ko, K2, K3, K, zo) 
+ (JV*(K4 - K, zi)iV*(K2 - Ko, z2))i^™,„,p,,(Ki, Ko, K3, K, zq) 
-(7V*(K4 - K, zi)iV(K3 - Ko, z2))i^,„,„.p,,(Ki, K2, Kq, K, zo) 
+ (iV*(K4 - K, zi)iV*(K - Ko, Z2))F,„,„,p.,(Ki, K2, K3, Kq, zo) 
d^Ko 

Xl^i;^dz2 dz. 



In Appendix |X] it is shown that 

/ / '(iV(Ki,Z2)A^*(K2,zi)) dz2 dzi =27r2dz5(Ki-K2)$i(Ki), 

J Zn J Zn 



(46) 



(47) 



where we set z — zo — dz. 

We now use Eq. (H71) to simphfy Eq. and then take the limit dz 0, to turn it into a differential equation 
again, 



(Ki,K2,K3,K4,z) 



2k 



(|Kip - |K2|2 + IK3P - |K4|2)F„,,„^p,5(Ki,K2,K3,K4,z) 



47r2 



2fc2F„,„,p,,(Ki,K2,K3,K4,z) j $i(K 

$i(K)F™,„,j,,,(Ki ~ K, K2 - K, K3, K4, z 
$i(K)F™,„,p,,(Ki, K2 - K, K3 - K, K4, z 



$i(K)F™,„,p,,(Ki - K, K2, K3, K4 - K, z 
j $i(K)F™,„,p,,(Ki, K2, K3 - K, K4 - K, z 
j $i(K)F„,„,p,,(K + Ki, K2, K3 - K, K4, z 
-e j $i(K)F™,„,p,,(Ki, K + K2, K3, K4 - K, z 



4^ 

d^K 

47r2 

d^K 

4^ 

d^K 

4^ 

d^K 

4^ 

d^K 

47r2 ■ 



(48) 



If one substitutes Eq. into the z-derivative of Eq. ([57]) . one would obtain a first order differential equation for 
the density operator. However, we are interested in the transformation of the individual density matrix elements. 
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B. Extraction of matrix elements 



To express the transformation of the density operator due to the infinitesimal propagation through a turbulent 
atmosphere in terms of the density matrix elements, we extract the matrix elements from the density operator using 
the trace 



dzPu,v,r,s{^') = tTi&ce {dzp{z)\v)\s){u\{r\} , 
where the operator that selects a particular matrix element in the 0AM basis is given by 



m{u\{r\ = \j |K8)|K6) F„%_,,,(K5,K6,K7,K8,z)(K7|(K5 



d^j^5 d^gfi d?K^ d^Ks 
Air^k.K Air^k^a Air^k^v Air^k^; 



(49) 



(50) 



The factor of a quarter in Eq. (|50p comes from the fact that, propagating through the same medium, the two photons 
are indistinguishable. As a result there are four different ways in which the states can be contracted on each other. 
This gives rise to a symmetry factor of 4, which implies that the same terms are counted several times. To remove 
this over counting, one needs to insert the factor of a quarter in Eq. (j50p . 

We now substitute the derivative of Eqs. (|37p with respect to z and ((50)) into Eq. (|49)) to obtain 



dzPu,v,r,s{z) = trace < ^ Pm,n,p,q / | Ki) IK3) 9^i^„,„^p^g (Ki , K2 , K3 , K4, z) 

m,n,p,q 



X(K2|(K4 



d'^Ki d^K2 d^Kg d^Ki 
47r2fc2i 47r2fc22 47r2fc23 Air'^k 



Air'^k, 



Pm,n,p,q j 52Fm,n,p,g (Ki , K2 , K3 , K4, z) 

d^Ki d^K2 d^i^g d^Ki 



m,n,p,q 

xF* (Ki,K2,K3,K4,z) 



47r2fc2i 47r2fc22 47r2fc23 An'^kzi ' 



(51) 



where the last expression is obtained because of the orthogonality of the momentum basis, given in Eq. (|24p . The 
factor of 4 that comes from the multiple ways in which the momentum states can be contracted on each other removed 
the factor of a quarter in the last expression. Substituting Eq. (|48l) into Eq. (|5T|) . we obtain 



^zPu^v^r.s (^) 



2k 



|Ki| 



IK2P 



IK: 12 



^ ^ Prn^n.p.q 
m,n,p,q 

-2fc2i^„,„,p,,(Ki,K2,K3,K4,z) J $i(K 

J <i>i(K)F„,„^p,,(Ki - K, K2 - K, K3, K4, z 
$i(K)i^„,,„^p^,(Ki, K2 - K, K3 - K, K4, z 
<i>i(K)F„,,„^p,,(Ki - K, K2, K3, K4 - K, z 
$i(K)i^„,,„^p^,(Ki, K2, K3 ~ K, K4 - K, z 
$i(K)^^„,,„^p,,(K + Ki, K2, K3 - K, K4, z 
-fc2 J $i(K)F™,„,p,,(Ki, K + K2, K3, K4 - K, z 
xi^* .,,,(Ki,K2,K3,K4,z) 



31 - |K4|')i^™,„,p,q(Ki,K2,K3,K4,z) 

d^K 



47r2 

4^ 
d^j^ 
47r2 

d2j^ 

47r2 
4^ 

d2j^ 

4^ 

d2^" 

4^_ 

d^i^i d2i^2 d^K3 d^K^ 



An'^kzi 47r2fc22 47r2/c23 47r2fc24 ' 



(52) 



The expressions in Eq. (j52p can be further simplified by using the orthogonality of the momentum space wave functions 
of the 0AM basis given in Eq. 
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C. Final expression 



After the final simplification of Eq. ([52)) we obtain the 
IPE as a set of first order differential equations given by 

(^zPu,v,r,s — Sm.uPm,v^r,s ^nPu^ri^r^s 
'^Sp^r Pu^v ^p^s Ss^qPu.v,r^q 
'^-^v,u,n.mPm.n,r.s ~t~ -£'D,r,n,pPu,n,p,s 
^"-^s,^^,q,mPm,^),r,Q ~t~ Ls^r.q,pPu,v^p^q 
Lm.r,u.pPm,v.p.s I^s,n^q^vPu,n^r,q 
— '2Pu,v,r,sLT, (53) 

where repeated indices imply summation. The quantities 
in Eq. ([53| are defined as follows 



Sx,y{z) 



2k 



|K|2g,(K,z)G*(K,z) 



Lr^k^ J *i(K) 



47r2 ' 



(54) 



(55) 



and 



im,n,n,.(^) = A:2y"$i(K)W,„,„(K,z)W^„%(K,^) 
with 



W^JK,z)= J G,(Ki,z)G;(Ki-K,z) 



d^Ki 



(57) 



The first four terms of Eq. ((53)) are the non-dissipative 
terms, representing the free-space propagation process. 
The last seven terms of Eq. ()5^ are the dissipative terms, 
representing the scattering of 0AM modes into other 
0AM modes due to scintillation caused by the turbu- 
lence. 



where p, q, and w are used to generated the spectrum 
of a particular LG mode, a and b are normalized spatial 
frequency components related to kx and ky via 



27ra 

/C^; ' 

do 



2TTb 

do 



and 



C2 dl ^ 



1/4 



(59) 



(60) 



with Wo and do being the centre frequency and beam 
radius, respectively. 

When the generating function in Eq. (|55)) is substituted 
into Eq. ((54)) the fc^-factors cancel, leaving the same ex- 
pression that was evaluated in Q . 

The integral in Eq. ((57| represents the correlation be- 
tween the momentum space wave functions of the 0AM 
modes. In this case the fc^-factors do not cancel, but 
leave a factor of 



1/2 



d2-A2(ai-a)2-A2(6i-fe)2 



dl 



1/4 



, (61) 



where A is the wavelength associated with ujq. In the 
paraxial limit (A <C do), this factor becomes 



1/2 



1 



(62) 



Therefore, one can neglect this factor, rendering the re- 
sulting expression in the same form as that which was 
evaluated in 

Hence, in the paraxial limit all the integrals in 
Eqs. ((M)) - ((37)) are the same as those that were con- 
sidered in 1^. The current formalism that starts from 
Lorentz invariant definitions of the quantum states does 
not produce expressions that are significantly different 
from those that were obtained in [^]. 



VI. SOLVING THE INTEGRALS 



VII. CONCLUSIONS 



The solution of the integrals in Eqs. ((54 )) - ((57)) have 
been considered in 01 . Two of these expressions, 
Eqs. ((M)) and ((57)) . are now slightly different due to the 
presence of the 1/fcz-factor. Moreover, to satisfy the 
orthogonality condition for the momentum space wave 
functions G„i(K,z), the generating function for the LG 
modes that was proposed in [1] , now needs to contain an 
additional factor of the square root of fc^. The expres- 
sion for the generating function of the spectra of the LG 
modes is therefore given by 



nk 



1/2 



■ exp 



«7r(a -|- ib)p in^a — ib)q 



7r2(a2 + 62)(i 



' w 
w - 



it- 



1 + w 
iwt) 



1 + w 



(58) 



An expression is obtained for the evolution of the den- 
sity matrix elements of a spatial mode entangled bi- 
photon propagating through the same turbulent atmo- 
sphere. The expression represents an infinite set of cou- 
pled first order differential equations, each containing 
non-dissipative terms associated with free-space propaga- 
tion and dissipative terms associated with modal scatter- 
ing due to the scintillation. Among the dissipative terms 
are included terms associated with the cross correlation 
between the two photons, coming from the fact that they 
are indistinguishable photons propagating through the 
same medium. 

The derivation starts from a manifestly Lorentz co- 
variant formulation of the quantum states and then fol- 
lows a number of steps that eventually explicitly break 
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the Lorentz invariance. This leads to a final expression 
that differs from the previous expression in Q in that 
some of the integrals (the phase space integrals) in the 
new expression contain factors of l/fc^. The only place 
where these factors could potentially have an effect is in 
the correlation between different momentum space wave 
functions for the different 0AM modes. However, in this 
case it is shown that the l/fc^-factors can be ignored in 
the paraxial limit. 

Hence, although the formal expression that is obtained 
from an initially Lorentz invariant formulation differs 
from what is obtained from the more traditional formula- 
tion, the effect of this difference on measurements is not 
expected to be significant or even observable. 



Appendix A: Ensemble average 

The evaluation of the ensemble averages of products of 
the random functions N{K,z) in Eq. pS)) . is discussed 
in an appendix in For the sake of convenience this 
appendix is included here, with only minor changes. 

As mentioned in Sec. IIV[ the refractive index fluctua- 
tion produced by a turbulent atmosphere is small com- 
pared to the average refractive index of air, 5n <^ 1, 
which leads to the fact that one can separate the propaga- 
tion through a turbulent atmosphere into two parts: free- 
space propagation and random phase modulations. The 
random phase functions for the latter step are obtained 
by integrating the refractive index fluctuation through a 
thin sheet of atmosphere along the propagation direction, 

/•zo+Sz/2 

0{x,y) ~ k Sn{x,y, z)dz 

Jzo-Sz/2 

« k 5z 5n{x,y, Zq), (A1) 

where, in the last line we took the limit 5z 0. Re- 
placing the refractive index fluctuation with its Fourier 
expansion, we obtain 

e{x,y,zo) = kSz / exp[-i{k^x + kyy + k^zo)] 



xA^(k) 



d^fc 
(27r)3 



(A2) 



where iV(k) is the three-dimensional spatial spectrum of 
the index fluctuation. We now define a two-dimensional 
spectrum for the accumulated index fluctuation A^(K, z) 
over a thin sheet of atmosphere 

/dk 
exp{-ih.z)Nik) (A3) 

which depends on the z position of the thin sheet. The 
three-dimensional spectrum of the refractive index fluc- 
tuation can be expressed in terms of its three-dimensional 
power spectral density, which follows from the autocor- 
relation function of the index fluctuation and which rep- 
resents the model for the turbulence, 



Nik) = x(k) 



$o(k) 



1/2 



(A4) 



where x(k) is a normally distributed random complex 
spectral function and is its spatial coherence length 
in the frequency domain. The latter is inversely propor- 
tional to the outer scale of the turbulence. Since the re- 
fractive index fluctuation Sn is an asymmetric real- valued 
function, we have that x*(k) = x{—k). Furthermore, the 
autocorrelation function of the random function is 



(X(ki)r(k2)) = (2^Afc)3j3(ki - k2). 



(A5) 



In Eq. (|46p we flnd ensemble averages inside double 
^-integrals. Substituting Eqs. (|A3|) and (|A4[) into such 
an ensemble average and using Eq. (jA5[) to evaluate the 
ensemble average, one obtains 

f f ' (7V(Ki, Z2)N*{K2, Zl)) dz2 dzi 
J Zn J Zn 



I Zq J Zf) 

\2 



(2^)^5(Ki - K2) 



X exp [ikz{zi — Z2)] dz2 dzi 



$o(ki) 

dfc. 
27r 



(A6) 



We set z = Z{) + dz and evaluate the two z-integrals 

(•ZQ+dz pzi 



exp [ikz{zi - Z2)\ dz2 dzi 



1 — coslk^dz) sm(kzdz) — kzdz , , 
+ « — • (A7) 



The power spectral density $o(ki) is always even in k^. 
Therefore, the imaginary part of Eq. (jA7|) . being odd in 
kz, does not contribute to the final expression, 

r r{N{K,,Z2)N*{K2,z,)) dz2 dzi 

J zo J za 

^ (2^)25(Ki - K2) j $o(ki) 



1 — cos{kzdz) 



1.2 



dfc. 
27r 



(A8) 



Due to the fact that the refractive index variations are 
very small, the light that propagates through the turbu- 
lent atmosphere remains unchanged over distances much 
longer than the correlation distance of the turbulent 
medium. One can therefore assume that dz is much 
larger than this correlation distance. As a result the 
function inside the square-brackets in Eq. (lASp acts like 
a Dirac delta function, so that one can substitute kz = 
in $0 E^nd pull it out of the fc^-integral. The integral can 
then be evaluated to give 

f ' / {N{Ki,Z2)N*{K2,z)) dz2 dz 

J Zq J Zq 

= 27r2dz,5(Ki -K2)$i(Ki), (A9) 
where we defined <l'i(Ki) — $o(Ki,0). 
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